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Abstract. We show that the emergence of criticality in the locally-defined Bak-Sneppen model corresponds 
to separation over a hierarchy of timescales. Near to the critical point the model obeys scaling relations, 
with exponents which we derive numerically for a one-dimensional system. We further describe how the 
model can be related to the glass model of Bouchaud [J. Phys. I France 2, 1705 (1992)], and we use this 
insight to comment on the usual assumption of stationarity in the Bak-Sneppen model. Finally, we propose 
a general definition of self-organised criticality which is in partial agreement with other recent definitions. 
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1 Introduction 

The concept of self- organised criticality or SOC was origi- 
nally devised by Bak, Tang and Wiesenfeld as an explana- 
tion for the apparent ubiquity of scale invariant systems in 
nature Ji10,0,EIh,o|. Working on the assumption that such 
systems are 'critical' in the sense of a continuous phase 
transition, Bak et al. proposed the existence of a class of 
non-equilibrium models that become critical purely under 
their own dynamics. This contrasts with other equilibrium 
and non-equilibrium systems, where at least one control 
parameter must be set to a particular value before the 
critical state is attained [0,||] , which seems unlikely in the 
absence of human intervention. Despite these early claims, 
it soon became clear that many SOC models do in fact re- 
quire parameter tuning, but they had been defined in such 
a way that the tuning had been carried out implicitly. For 
instance, the sandpile model, which is the canonical SOC 
system, is usually described in a way that implicitly as- 
sumes an infinitesimal driving rate; for finite driving, the 
model is no longer critical P,pQ] . The admission of implicit 
parameter tuning has allowed the sandpile model to be re- 
lated to more conventional non-equilibrium systems, and 
for analytical techniques previously established in other 
fields to be employed (see e.g. [ p] , [l^ ,[l^) 

Recently, attention has focused on implicit parameter 
tuning in a subclass of SOC systems known as extremal dy- 
namical systems, so-called because they are driven at the 
location of the minimum (or maximum) of some spatially 
varying quantity |l4,n5|. It has already been observed in- 
dependently by Sneppen |l^ and Vergeles ||l^ that the 
Bak-Sneppen model, which is the simplest and perhaps 
best understood extremal dynamical system [p^,p^,[l9|, 
implicitly assumes that a single temperature-like param- 



eter is set arbitrarily close to zero. However, in our view 
the mechanism underlying this process was never clearly 
identified, nor was the inevitability of the parameter tun- 
ing properly stressed. 

The aims of this paper are fourfold. Firstly, we attempt 
to provide as clear an explanation as possible that implicit 
parameter tuning in the Bak-Sneppen model, rather than 
just being possible, is in fact inevitable if one is to have 
a well-defined physical system. By 'well-defined' we mean 
that each element evolves according to the state of only a 
finite number of other elements, as opposed to the global 
driving rule in the model's original invocation. This relates 
to recent work on the sandpile model M, but was never 
properly discussed in the previous work on this model [16l 
p7[ . Note that throughout this paper we follow [^ in refer- 
ing to the setting of a parameter to zero as 'fine tuning,' 
although it could be argued that this is just an absolute 
separation of temperature scales rather than genuine fine 
tuning. Secondly, we identify the fundamental mechanism 
underlying this process to be a hierarchy of timescales 
that diverge relative to each other, which is similar but 
qualitatively different to the separation of timescales re- 
quired in other SOC models pO|. Thirdly, we show that the 
model can be related to a phenomenological glass model 
of Bouchaud |gl| , and argue that the Bak-Sneppen model 
does not reach a statistical steady state in low dimensions, 
in contrary to the widely held belief that it reaches sta- 
tionarity after an 'extensive transient' p8| . Finally, we use 
the insight gained from our work to propose a new way of 
categorising self-organised criticality, and comment on its 
relationship with other recent definitions. 

This paper is arranged as follows. In section 2 we care- 
fully reconstruct the Bak-Sneppen model from first prin- 
ciples, and show how a properly defined model with only 
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local interaction rules demands implicit parameter tuning, 
as in the sandpile model |^. In section 3 we show that at 
finite temperatures the model can be related to a simple 
glass model, and argue that even at its critical point, the 
model does not reach a steady state. A mean field descrip- 
tion similar to that adopted by Bouchaud is described and 
solved in section 4, and we comment on the true meaning 
of self-organised criticality section 5. 



2 Reconstruction of the Bak-Sneppen model 

We begin by reconstructing the Bak-Sneppen model from 
first principles, following the original arguments in |18| . 
Although some aspects of this have already been touched 
upon p^,pT[, the emphasis here is to demonstrate utterly 
unambiguously that the implicit parameter tuning is un- 
avoidable, rather than just possible. 

The local Bak-Sneppen model is defined as follows. 
The system consists of N elements, each of which is as- 
signed a barrier Ei , i = 1 . . . N . In the model's original bi- 
ological context the Ei represent fitness barriers, but they 
could just as easily correspond to energy barriers on a po- 
tential energy landscape, for example. The values of the 
barriers are drawn from the time-independent prior distri- 
bution p{E), which is assumed to have no delta- function 
peaks so that there is a vanishing probability of two dif- 
ferent elements having the same value of E. The system 
evolves according to two rules. Firstly, each element be- 
comes activated at a rate 



-Ei/T 



(1) 



where the constant parameter T > has the same units 
as E. An activated site i is assigned a new barrier Ei drawn 
from p{E), corresponding to a shift to a new metastable 
state with a new barrier height. Note that the activation 
rate for any given element is independent of the state of the 
rest of the system, so that this activation rule is strictly 
local. 

Secondly, for every activated element, another z are 
chosen and also assigned new barrier values. This interac- 
tion term mimics some form of strong coupling between 
the elements, in the sense that one element changing state 
drastically alters the barriers of z other elements. The 
way in which the z interacting elements are chosen de- 
pends upon the spatial structure of the system. If the el- 
ements are arranged on a regular lattice, then the most 
common rule is to update the barriers of all of the near- 
est neighbours of the activated element. Hence z is just 
the lattice coordination number. Alternatively, the z in- 
teracting elements may be chosen at random from the re- 
maining A^ — 1, with the connections between the elements 
being randomised anew for every activation event. In this 
case the system has no spatial structure and the symbol 
K has often been used, where K = z + 1 ]2^,p3[. 

The system behaviour changes qualitatively as the sin- 
gle parameter T > is varied. These different regimes are 
discussed in turn below. 



T -^ 0^; In the limit of infinitesimal T, the activation 
rates e^^'' "^ for different Ei diverge relative to each other. 
That is, the element with the smallest barrier will become 
active on one timescale, the one with the second smallest 
barrier becomes active on another, much longer timescale, 
and so on. Thus with probability one the first element to 
become active will be that with the smallest barrier (which 
is always unique for a finite set of non-degenerate {Ei}). 
This is the way in which the Bak-Sneppen model is usually 
defined; indeed, this "exponential separation of timescales" 
[18| was originally used to justify the extremal dynamics. 
Although separated timescales arise in other SOC mod- 
els [^0| , to our knowledge this is the first time a hierarchy 
of timescales has been explicitly demonstrated. 

T small but finite: The strict separation of timescales is 
lost for finite T and every element has a non-vanishing 
probability of being the first to become active, so the dy- 
namics are no longer extremal. We now demonstrate that 
the model is not critical in this regime. This claim is sup- 
ported by the results of numerical simulations described 
below, which in turn are supported by mean field analysis. 
Let Pi denote the probability that element i becomes 
active before any other element in the system. Since the 
time until activation follows an exponential distribution 
with mean e^'/"^, it is straightforward to show that 



-Ei/T 



EU e-^^/^ 



(2) 



^N 



which obeys J2i=iPi = 1- I^ this notation, the T ^ 0+ 
limit corresponds to pi<- — » 1 for the element i* with the 
smallest barrier, and pj -^ for all j ^ i*. However, < 
Pi < 1 for all i when T is finite. 

The algorithm employed in the simulations was as fol- 
lows. N elements were placed on a one dimensional lattice 
with periodic boundary conditions. Each element was ini- 
tially assigned a barrier drawn from the uniform prior dis- 
tribution p{E) = {1 for < i? < 1, otherwise}, although 
we expect the same qualitative behaviour for other p{E). 
For every iteration step, a single element was made ac- 
tive according to the probabilities pi given in (g). The 
active element and both of its nearest neighbours were 
then assigned new barriers, so the number of interacting 
elements z = 2 here. The timescale t was normalised to N 
activations per unit i, which differs from the usual Bak- 
Sneppen timescale only by the constant factor N . Distri- 
butions were not measured until the mean barrier height 
{E) = j^ X]i=i ^i appeared to reach a steady value when 
plotted against logj^p ^- Note that this is not a rigorous 
criterion for convergence and we cannot rule out the pos- 
sibility that long-range spatial correlations may still be 
growing. 

To decide whether or not the system is critical for any 
given value of T, we employed the usual method of extract- 
ing the spatial and temporal correlations from the simula- 
tions and checking to see if their tails are consistent with a 
power law fit. Apart from complications arising due to fi- 
nite size effects, anything other than power law behaviour 
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Fig. 1. Log-log plot of the distribution of the .jump size \x\ 
between successive active sites for an N = 10'* ring. From bot- 
tom to top on the right hand side, the lines refer to T — > 0^ 
(thick line), 0.001, 0.002, 0.005, 0.01, 0.02 and 0.05 (thin lines), 
respectively. The thick line has a slope of —3.20 ± 0.05. (In- 
set) The mean barrier {E} against logj^Q t for T — 0'^ (thick 
line), 0.01, 0.02, 0.03, 0.04 and 0.05 (thin lines, top to bot- 
tom). Each unit t corresponds to TV activations. Data collection 
commenced at logj^g i = 4 in this and all subsequent figures. 



signifies a characteristic scale and a non-critical system. A 
convenient measure of spatial correlations for this model 
is the number of lattice sites between two successive ac- 
tive elements. The distribution of these 'jump sizes' jxl for 
different T are plotted in Fig. 0. For T ^ 0+ we find that 
-fjump(|a;|) ~ la^I"'^ with tt — 3.20 ± 0.05, in accord with 
the known value 3.23 ±0.02 [Q. However, Pjump(|a^|) lev- 
els out at a constant T-dependent value Pjump(oo) ~ ^(T) 
for finite T, indicating that the system is not critical. As 
demonstrated in Fig. | for an iV = lO'^ system, A{T) ~ T" 
with a = 3.0 ± 0.2, showing that the critical point is in- 
deed at T — > 0+. This is consistent with the results of 
Sneppen p6[ , but our data appears to be much smoother, 
allowing for a more precise evaluation of the exponent. 

The temporal correlations are quantified by Pret (i) , the 
distribution of times t since the currently active site was 
last active, which is plotted in Fig. 0. Again the data for 
finite T is clearly not power law, and furthermore the 
data for small T obeys a scaling function of the form 
Pvetit) ~ t-°'%l^{TPt) with ip{y) -^ (const) as y ^ 0, with 
the exponents a — 1.58 ± 0.01 and /3 = 3 ± 0.1, as demon- 
strated in the inset of Fig. 13. 

We interpret the loss of criticality for finite T as fol- 
lows. For T — > 0+ the location of the active site jumps 
around the system in the highly correlated manner char- 
acteristic of the critical state. By contrast, when T is finite 
there is a non-zero probability that elements arbitrarily far 
from the active site will become activated at the next time 
step. Thus the active site can make large jumps to uncor- 
related regions of the system, which we relate to the loss 
of criticality. 




Fig. 2. Log-log plot of the limiting value lim|j;j^oo Pjump(|a;|) 
as a function of T for A'^ — lO**. The dashed line has a slope 
of 3. 
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Fig. 3. Log- log plot of Prct(t), the distribution of times since 
the currently active site was last active, for different T on an 
A'^ = 10* ring. From bottom to top in the middle of the graph, 
the lines refer to T ^ 0+ (thick line), 0.001, 0.002, 0.005, 0.01, 
0.02 and 0.05 (thin lines), respectively. The thick line has a 
slope of -1.58±0.02. (Inset) Scahng plot of i^'^**Prot(t) against 
T^t for T = 0.001 (solid line), 0.002 (dotted line) and 0.005 
(dashed line). 



T = 0{E): Although there is no longer any question of 
criticality far from T — 0"^, it is possible to relate the 
model to a glass model for finite T. This will be discussed 
fully in section 3. 

T — > oo.' In this limit, every element has the same activa- 
tion probability pi = 1/iV independent of E[ , so P{E, t) = 
p{E) for all t. 
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3 Mean field analysis 

In this section we show how the results from the one di- 
mensional simulations are in qualitative agreement with 
the solution to the mean field model. The mean field model 
described here is essentially that of Bouchaud et al. with 
an extra interaction term pq | . 

Let P{E, t)dE be the proportion of barriers in the 
range [E, E + dE) . By adopting the usual mean field sim- 
plification of random nearest neighbours, it is then straight- 
forward to show that P{E, t) evolves in time according 
to M 



dPjE, t) 
dt 



Jo 



r, P{E, t) 



-zP{E,t) + {z + l)p{E) 



(3) 



This equation can be justified by noting that P{E, t) de- 
creases when an element changes its barrier value, which 
occurs either when it becomes active, or when it is selected 
as one of the z interacting elements. These two processes 
are described by the first and second terms on the right 
hand side of (||), respectively, where the prefactor to the 
first term is just the continuum analogue of (||). Conser- 
vation of probability is ensured by the third term, which 
corresponds to the z -\- \ new barriers drawn frorn p(E). 
We note that the treatment of Veregles is similar [0 but 
with a poorly defined timescale, resulting in factors of N 
remaining even after taking the continuum limit. 

We have solved ^ for the uniform p{E) in the limit 
i — > oo. The expression for general T is not very instruc- 
tive, but for small T it simplifies to 



P{E,oo) 



1 



1^^-{E-E^)/T 



(4) 



with Ec^ -^ 



As T ^ O'^ the exponential in (^ either 
blows up or decays depending on whether E is less than 
or greater than E^ , respectively. Thus P{E, cxo) converges 
to the step function ^^ 6{E — Ec), in accord with the 
known solution of the mean field Bak-Sneppen model |0, 




0.4 0.6 

Barrier value E 



Fig. 4. The barrier distribution P{E) for the same runs as in 
Fig. hi From top to bottom on the right hand side, the lines 
refer to T = 0+ (thick line), 0.01, 0.02, 0.03, 0.04 and 0.05 
(thin lines) , respectively. (Inset) The corresponding mean field 
predictions from (M, with 2 = 2. 



dt_ 
dr 



p-E'/T 



P{E\ t) dE' 



(5) 



The mean field equation (H) then becomes 



dP{E, t) 



/•OO 

,-E/Tp^^^ r)+p{E) / c-^'/^P{E', t) dE', 
Jo 

M 



,|24|. However, there is no such discontinuity for finite 
T and P{E, oo) is smoothly varying for all < i? < 1, in 
qualitative agreement with the simulation results in Fig. ^ 
(although note that Ec is larger in the one dimensional 
case). Note that ^ can also exhibit glass-like behaviour 
for T < 1, as fully described in |2q|. 



4 Glassiness and the assumption of 
stationarity 

Although the model studied here was described as the 
Bak-Sneppen model extended to finite temperatures, it 
can equally be viewed as the simple glass model of Bouchaud 
with an extra interaction term pi| , |2q ,pq,p7[. The mapping 
to Bouchaud's model is acheived by the two stage process 
of first 'switching off' the interactions, i.e. setting z — 0, 
and then mapping to a timescale t which obeys 



which is the master equation to Bouchaud's model [E7|. 
This alternative interpretation becomes particularly rel- 
evant for va]ues of T comparable to the expected bar- 
rier height E = jEp{E) dE, as it is in this regime that 
Bouchaud's model predicts a glass transition for a certain 
class oi p{E). 

Since we have already described the relationship be- 
tween these two models in some detail elsewhere [ p6[ , we 
will not repeat those results here. However, it highlights 
what may be a very serious problem common to almost all 
previous treatments of the Bak-Sneppen model. In short, 
we believe there is already sufficient evidence that the 
Bak-Sneppen model does not reach a statistical steady 
state, just as Bouchaud's model does not reach stationar- 
ity for low temperatures. Clearly this is contrary the the 
widely held assumption that it does, so this point merits 
further discussion. 

It is widely known from glass theory that, to truly test 
a system's stationarity, it is not sufficient to measure func- 
tions of only one time variable. Instead one must measure 
a two-time correlation function C(t+tvj, tw), which is some 
measure of correlation between the state of the system at 
times tw and iw + 1, and show that it does not depend on 
the waiting time iw for large iw |^, ^ ■ Within the current 
context, almost all of the data commonly recovered from 
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simulations of the Bak-Sneppen model (critical exponents, 
P{E,t), etc.) are all functions of at most one time vari- 
able. For example, simply observing that P{E,t) appears 
to approach a limiting distribution Poo (E) does not show 
that the system has become stationary. Nonetheless it ap- 
pears that criteria similar to this are usually employed to 
check for a statistical steady state. 

To our knowledge, there has only been one instance 
when two-time correlations have been actively searched for 
in the Bak-Sneppen model, and this is the numerical work 
of Boettcher and Paczuski in one and two dimensions pO| . 
Remarkably, they found that the system behaviour clearly 
depends on t^ for all t^ they measured. This situation is 
commonly referred to as aging and is a clear indicator 
of non-stationarity (ie. loss of time translational invari- 
ance). That aging implies non-stationarity is trivial; the 
only question that remains is, are Boettcher et aL's re- 
sults asymptotic, or is stationarity recovered at some very 
late time r which is beyond attainable simulation times? 
Clearly this can never be answered by numerical simula- 
tions alone and some form of analytical treatment would 
be desirable. However, we have shown elsewhere that the 
mean field model does not exhibit aging |26[ , so analysis 
would have to limited to the difficult case of finite dimen- 
sional systems. This is an important issue whose resolution 
may help guide attempts to find an exact solution to the 
model, and further work would be desirable. 



scale, i.e. they are absolute rather than relative points. If 
any one of these points also happens to be a critical, then 
it is conceivable that the model could be placed at its crit- 
ical point 'by accident' and thus be erroneously referred 
to as 'self-organised' critical. It is becoming increasingly 
clear that this is precisely the case is all models currently 
referred to as SOC. 

In summary, we have elucidated the mechanism behind 
the approach to the critical point in the local Bak-Sneppen 
model. The underlying feature is the existence of a hierar- 
chy of timescales that become separated as a single con- 
trol parameter approaches zero. We have used the insight 
gained to propose a new definition of SOC that encom- 
passes all cases of implicit parameter tuning observed so 
far. Furthermore we have suggested that the Bak-Sneppen 
model may be non-stationary. This is based on the rela- 
tionship with a glass model for finite T, and suggests that 
other SOC models may also exhibit interesting behaviour 
far from their critical points. We welcome study of these 
and related questions. 
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5 Discussion and summary 

It has been suggested in g that SOC systems correspond 
to absorbing state phase transitions reached in the limit 
of infinitesimal driving. Although this is almost certainly 
true, in our view this only explains why SOC systems are 
critical, not why they are self-organised. That is, it does 
not address, in sufficiently general terms, how a model 
can be placed at a particular point on its phase diagram 
without explicit parameter tuning. Note that this is a sep- 
arate issue as to whether or not the point also happens 
to be a critical point. An alternative definition of SOC 
has addressed the self-organising process on more general 
terms [^,32|, but gives httle insight into what class of 
systems need to be critical to reach a statistical steady 
state. 

In light of the work presented in this paper, we now 
propose a broad definition of SOC that is relevant to the 
phase diagram approach adopted in H and throughout 
this paper. We suggest that SOC corresponds to that class 
of models that have a critical point at a privileged point on 
their phase diagrams. By 'privileged' we mean any point 
at which every parameter takes a value that has some spe- 
cial physical significance. For instance, a positive definite 
parameter such as the rate of driving or temperature has 
two privileged values, 0"*" and oo, and indeed all implicit 
tuning identified thus far do seem to require infinitesimal 
driving or temperature |lq,EG|. Similarly a conservation 
parameter has special points corresponding to 100% con- 
servation and 100 % dissipation. Note that such points are 
scale invariant in that they do not depend upon the chosen 
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